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Integral equation theory is applied to a coarse-grained model of water to study potential of mean 
force between hydrophobic solutes. Theory is shown to be in good agreement with the available 
simulation data for methane-methane and fullerene-fullerene potential of mean force in water; the 
potential of mean force is also decomposed into its entropic and enthalpic contributions. Mode cou- 
pling theory is employed to compute self-diffusion coefficient of water, as well as diffusion coefficient 
of a dilute hydrophobic solute; good agreement with molecular dynamics simulation results is found. 
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I. INTRODUCTION 



Hydrophobic interactions play an important role in a wide variety of phenomena, including collapse of hydrophobic 
polymers in water, micelle formation, and protein folding. 0, Q As a result, this problem has been actively studied 
theoretically. The interactions between hydrophobic solutes in water have been investigated using both computer 

simulations and mean field theoretical methods. The former approach produces exact results for a given microscopic 
model, but at a significant computational cost, especially for large-size solutes, which require simulating a very large 
number of water molecules. The latter approach is much less computationally demanding, but inevitably involves 
approximations. 

One possible way to reduce computational cost associated with simulation studies is to employ simplified coarse- 
grained solvent models. [IH 0| Of particular relevance to the present study is the recent work of Shell and 
co-workers, [23;, 24] who have shown that several thermodynamic and dynamic anomalies of pure water, as well as 
various features of hydrophobic interactions can be reasonably well reproduced with a coarse-grained water model 
based on an isotropic pairwise additive "Lennard- Jones plus Gaussian" (LJG) interaction potential. While the work 
of Shell and co-workers employed simulation-based techniques, their microscopic model based on an isotropic pair 
potential is ideally suited for approaches such as i nteg ral equation theory and mode-coupling theory (MCT), which 
would reduce computational expenses even further. (25l - l27| 

Shell and co-workers have found that in order to properly reproduce the behavior of potential of mean force 
(PMF) between two hydrophobic solutes in water (e.g. its variation with the solvent temperature), the parameters of 
the coarse-grained LJG potential must be taken to be dependent on the thermodynamic state point of the solvent. 
Nevertheless, once the potential is parameterized for a given density and temperature of water, it can be employed to 
study the hydrophobic interactions between various solutes, and here is precisely where the computational efficiency 
of the integral equation approach could be very helpful. However, as mentioned above, the theory is necessarily 
approximate and its accuracy needs to be tested against the simulation data. The goal of the present study is to 
perform a comparison between theory and simulation in obtaining structural and dynamical properties of both pure 
water and dilute solutions of hydrophobic solutes. 

The remainder of the paper is organized as follows. In Section[H]we present our microscopic model and theoretical 
methods for calculating structural, thermodynamic, and dynamical properties. The results are presented in Section lHTl 
Section UVl concludes the paper. 



II. MICROSCOPIC MODEL AND THEORY 



A. Microscopic Model 

As a coarse-grained model for water, we consider a system comprised of spherical particles interacting via an 
isotropic LJG pair potential of the following form: [2^, [24| 



B exp 



(r - r ) 5 



(1) 



where e and a are the usual LJ parameters setting the energy and length scale of the interaction, while parameter B 
sets the strength of the Gaussian term, whose center and width are determined by ro and A, respectively. 

In studying hydrophobic interactions, we consider dilute solutions, where solutes interact with solvent particles 
via isotropic pair potential <f>(r) and with each other via another isotropic pair potential </> uu (r), particular forms of 
these potentials depend on the solute and will be specified for each solute in Section HTT1 



B. Structural Properties 

We employ integral equation theory 28] to compute the solvent-solvent (g s (r)) and the solute-solvent (g(r)) pair 
distribution functions and the solute-solute PMF (W(r)). 

The starting point of the calculation is the Ornstein-Zernike relation between the solvent-solvent total (h s (r) = 
g s (r) — 1) and direct (c s (r)) pair correlation functions: [28| 

h s (r) = c s (r) +p f dr'h s (r')c s (\r- r'\), (2) 
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where p is the bulk solvent number density. In order to solve the above equation, one needs an additional closure rela- 
tion between h s (r) and c s (r). Here we employ thermodynamically self-consistent hybrid mean spherical approximation 
(HMSA) closure: [I| 



,(r) = exp{-/3^ s (r)} 



exp[/(r)(fe s (r)-c s (r))]-l 
/(r) 



- h s (n) + c s (n) - 1, (3) 



where f3 — l/fc^T, and f(r) = 1 — exp(— ar). Note that this closure interpolates between the soft-core mean spherical 
approximation 30] (a — > 0) and the (hypernetted chain) HNC closure (a — s> oo).[29] 

The value of the parameter a is determined by a thermodynamic consistency condition that is based on equating 
the isothermal compressibilities, Xt ano - Xt> that follow respectively from the virial and compressibility routes: p9l ] 

<>«*>-- (3F 

{pkT X c T r l =1-9 I dfc s (r), (4) 



where the virial pressure P v is given by 

PP V =P~^f j ' dfr4>' s (r)g s (r), (5) 

with prime denoting differentiation with respect to the argument. 

For an infinitely dilute solute in a solvent, the solute-solvent total (h(r) = g(r) — 1) and direct (c(r)) pair 
correlation functions are also related via Ornstein-Zernike equation: 

h(r) = c(r) +pj df'h{r')c s (\f- f'\), (6) 

We employ the HNC closure to solve for h(r) and c(r): 

h{r) = exp[-^(r) + h{r) - c(r)] (7) 

The accuracy of this approach will be tested in the next Section by comparing theoretical results with simulations. 
Finally, for two dilute solutes in a solvent, the solute-solute PMF is composed of two terms: 

W{R) = <t> uu {R)+AW{R), (8) 

where the first term is the bare solute-solute interaction potential, while the second term is the solvent-mediated PMF 
given by the following exact relations: [3l| 



AW(R) = I F(R')dR', (9) 
where the excess mean force, F(R), is given by: 

F(R) = - I drV(j>{r)p{r;R). (10) 



In the above, p(f; R) is the conditional probability of finding the solvent particle at r given that one solute is at 
the origin, and the other solute is located at R. We compute this conditional probability from the anisotropic HNC 
closure, which has been thoroughly tested against computer simulations in our earlier work: (3ll . l32j 



p{r;R) = pexp 



p(4>{r)+4>Qf-8\j) + f df'c(\f-f'\)(p(r';R)-p) 



(11) 



From the solute-solute PMF one can easily obtain the solute-solute radial distribution function given by: 

9uu (R)=eM-l3W(R)). (12) 

As will be seen from the results presented below, for a typical hydrophobic solute such as methane, the solute-solute 
PMF typically displays two pronounced minima corresponding to the contact pair and the solvent-separated pair of 
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the solutes. @, [HI Concomitantly, the solute-solute radial distribution function is characterized by two corresponding 
maxima, with the second maximum (corresponding to the solvent-separated pair) bracketed by two minima located 
at Ri and i?2- In what follows, we will study the dependence of the methane- methane PMF and g uu (R) on the water 
density and temperature. In order to characterize the behavior of these functions by a single parameter, we will 
compute the equilibrium constant for the conversion of the contact pair into the solvent-separated pair:[7[ 

it 9uu(r)r 2 dr 
J Rl guu(r)r z dr 



C. Dynamical Properties 



We employ the MCT-based approach to compute the self-diffusion coefficient of the pure water, as well as the 
diffusion coefficient of a dilute solute. In this approach, the self-diffusion coefficient is obtained from the total friction 



D 



k B T 



(14) 



where m s is the mass of the water molecule, and the MCT expression for £ is comprised of binary and collective 
density contributions: [3r| 



C = Cb + C P 



(15) 



In what follows, we will be calculating self-diffusion coefficient at relatively high water densities, and therefore we 
neglect the hydrodynamic contribution to friction which arises from the coupling of the tagged solvent particle motion 
to the collective transverse current mode. [36] 

The binary term is given by the total time integral of the fast decaying binary component of the time dependent 
friction, which we model as a Gaussian function, whose parameters are chosen to reproduce the exact short-time 
behavior of the total time-dependent friction, ((t), up to the term of order t 2 :[33l| 



C b (t) =C(0)exp 



cm 2 

2C(0) 



with the initial-time value given by: 



C(o) 



4-7T/9 

3m, 



drr 2 g s (r)y 2 cf> s (r). 



(16) 



(17) 



The zero-time value of the second time derivative, ^(0) , can be similarly obtained from the fluid intermolecular 
potential and radial distribution function. [33l . |35| 

The density term arises from the coupling of the tagged solvent particle motion to the collective density mode: [36[ 



( P p(t)dt, 



with 



( PP (t) 



k B Tp 
6ir 2 m s 



dkk 4 c s (k) 2 [F s (k,t)F(k,t) - F°(k,t)F°(k,t)] 



(18) 



(19) 



where F{k 1 t) is the solvent dynamic structure factor, F°(k,t) = exp(— kp>Tk 2 t 2 /2m s ), and F s (k,t) is the solvent 
self-dynamic structure factor, for which we have adopted a simple Gaussian model: (33l [36| 



F s (k,t) = exp 



(20) 



We obtain the solvent dynamic structure factor from the continued fraction representation of its Laplace transform 
truncated at the second order: 1331 



F(k,z) 



z + 



z+- 



m 

Si(k) 
S 2 (k) 



(21) 
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where Si (k) is the initial time value of the i order memory function (MF) of F(k, t) . For the parameter r 1 (k) we use 
the expression due to Lovesey:[37j T _1 (fc) = 2y^(fc)/7r. The quantities Si(k) and 82(h) can be easily calculated from 
the first three short-time expansion coefficients of F(h,t); the microscopic expressions for the latter are well-known 
and will not be reproduced here. [H, [H, HI] 

Given that the self-dynamic structure factor in Eq. (I2TJ)) is a function of £, which, in turn, depends on F s (h, t) via 
Eq. (fl8)) and (|T9|) . the above set of MCT equations for £ needs to be solved iteratively and self-consistently. One could 
use a more accurate model for F 8 (k, t) in terms of the velocity time correlation function of a tagged fluid particle. [34j 
but our numerical calculations have shown that this does not change the results for D in a noticeable way. 

The diffusion coefficient of a dilute solute is obtained in the same way as the solvent self-diffusion coefficient, 
with the understanding that in the above set of MCT equations (f> s (r), g s (f), and c s (r) are replaced with <p(r), g(r), 
and c(r), respectively. Additionally, F s (h,t) now has the meaning of the solute self-dynamic structure factor, and 
Fg(h,t) — exp(— ksThH 2 /2m) is its inertial component, where m is the solute mass. Finally, m s is replaced with m 
in Eqs. 0, G3>, (X2)> and 420}. 

In the next section, we compare our theoretical results for structural and dynamical properties of pure water and 
dilute solutions of hydrophobic solutes with MD data and analyze the dependence of various computed properties on 
the solvent thermodynamic conditions. 



III. RESULTS 



A. Pure Water 



We start by presenting our results for the structural and dynamical properties of pure water. Parameters entering 
the LJG potential given by Eq. ([T]) are dependent on density and temperature and are taken from Ref. I23I where 
they were obtained from the relative entropy optimization method. In Fig. 1 we compare simulation results for the 
oxygen-oxygen radial distribution function of all-atom SPC/E water. [20j and the HMSA result for g s (r) of coarse- 
grained water described by the LJG potential. The two distribution functions are generally in good agreement except 
some discrepancies in the range between 3 and 6 A. The same discrepancies were seen in the earlier comparison of 
simulated g s (r) for SPC and LJG water, i.e. if one were to compare HMSA theory with simulations done for the 
coarse-grained LJG water (not shown), the two distribution functions would be nearly indistinguishable. 

The same level of agreement between theoretical and simulation results for the structural properties of coarse- 
grained water is obtained at other thermodynamic conditions, as can be seen from Fig. 2a, where we present MD[23| 
(symbols) and HMSA (lines) pressure isochores for two values of water density: p=0.9 g/cm 3 and 1 g/cm 3 , with 
theoretical results obtained from Eq. ([5]). One sees that the HMSA results for pressure are in excellent agreement 
with the MD data for both isochores throughout the entire temperature range studied, thereby confirming the accuracy 
of the theory in providing structural information for the coarse-grained water model. In particular, theory successfully 
reproduces the characteristic anomalous water behavior, with pressure decreasing as a function of T along the isochorc 
p=0.9 g/cm 3 at the lower end of the temperature range studied, subsequently passing through a minimum, and then 
growing with temperature. 

Having ascertained the accuracy of the HMSA theory in calculating structural and thermodynamic properties 
of pure water, we next use the structural data as input to compute water self-diffusion coefficient from the MCT 
approach. In Fig. 2b we present simulation [23| (symbols) and theoretical (lines) results for self-diffusion coefficient 
of the coarse-grained LJG water as a function of density along two isotherms: T=260 K and T=300 K. At the 
higher temperature, the self-diffusion coefficient decreases monotonically with increasing density, while at the lower 
temperature, anomalous behavior is observed, where the self-diffusion coefficient passes through a minimum around 
p=1.0 g/cm 3 , increases with density between 1.0 and 1.1 g/cm 3 , passes through a maximum, and then decreases with 
density. One sees that theory is in good agreement with simulation for both isotherms, and in particular captures the 
characteristic water diffusion anomaly well. 

Having considered structural, thermodynamic, and dynamic properties of pure coarse-grained LJG water, we next 
turn to dilute solutions of hydrophobic solutes in LJG water. 
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FIG. 1. Radial distribution function of pure water at T=300 K and p=l g/cm 3 . Symbols are the simulation results for the 
oxygen-oxygen radial distribution function of all-atom SPC/E water. [201] and the line is from HMSA theory for coarse-grained 
water described by the LJG potential. 



B. Methane Solutes in Water 



For the purpose of modeling dilute solution of methane in water, we have employed LJ potentials for both 
solute-solvent and solute-solute interactions: 



(r) = 4 £ll 



(22) 



and 



4>uu{r) = 4e„ 



(23) 



In order to compare our theoretical results with the earlier simulation study of methane solution in LJG water. [20L 1231] 
we have chosen the same values of e and a parameters as in the simulations: er us =3.45 A, e us =0.89 kJ/mol, a uu =3.73 A, 
and e u «=1.23 kJ/mol. 

In Fig. 3 we compare simulation results for the methane- water radial distribution function with all-atom SPC/E 
water, [2y and the HNC result for g(r) with coarse-grained LJG water. The two distribution functions are generally 
in good agreement except some discrepancies in the range between 4 and 7 A. As in the case of pure water, the 
agreement would be significantly better if one were to compare theory with simulation performed for methane in 
coarse-grained LJG water. This is confirmed by comparing theoretical and simulation results for the PMF between 
two dilute methane solutes in LJG water, which are presented in Fig. 4a for thermodynamic state point T=300 K and 
p=l g/cm 3 . One sees that theory is in excellent agreement with simulation throughout the entire range of separations 
between the two methane solutes, indicating that it accurately reproduces the anisotropic density distribution of LJG 
water molecules around two methane solutes at all separations considered. In particular, theory properly reproduces 




FIG. 2. (a)Pressure isochores for the coarse-grained LJG water model; symbols are from MD simulations,]^ and lines are 
from the HMSA integral equation theory, (b) Self-diffusion coefficient of the coarse-grained LJG water as a function of density 
along two isotherms; symbols are from MD simulations. [23l] and lines are from the MCT approach. 



the locations and the depths of the two minima in the PMF - the first one corresponding to the contact pair (located 
around 3.9 A), and the second one corresponding to the solvent-separated pair (located around 6.7 A); the contact 
pair minimum is deeper compared to the solvent-separated one. 

The agreement is equally good at the other two temperatures for which simulations were performed: T=280 K 
and T=320 K (not shown). On the basis of this information, we have decomposed the solvent-mediated PMF given 
by the integral equation theory into enthalpic and entropic contributions: |6l 

AS(T) ^W(T + ±n-AWiT-AT ) t (24) 

and 

AH = AW + TAS, (25) 

where AT=20 K in our calculations. The corresponding results are presented in Fig. 4b. In agreement with previ- 
ous simulation results, [ll| one sees that the solvent-mediated PMF at the location of the contact pair minimum is 
completely dominated by the entropic contribution, with enthalpic term very close to zero but slightly unfavorable 
(positive). At larger separations, the enthalpic contribution becomes favorable and the entropic term unfavorable; 
the solvent-mediated PMF at the location of the solvent-separated pair minimum is dominated by the enthalpic 
contribution. [ll| 

In order to present our results for methane-methane PMF at other densities and temperatures in a compact form, 
we have computed the equilibrium constant for the conversion of the contact pair into the solvent-separated pair from 
Eq. (|T3"]) . Theoretical results for K eq along two isotherms are presented in Fig. 5 together with the only simulation 
data point available[23] at T=300 K and p=l g/cm 3 . One sees that in agreement with earlier simulation work,Q the 
solvent-separated pair is more stable in terms of K eq at all densities and temperatures studied. For both isotherms, 
the equilibrium constant first increases with density, passes through a maximum between 1.05 g/cm 3 and 1.1 g/cm 3 , 
and then decreases with density. Except for the highest density considered (p=1.2 g/cm 3 ), the equilibrium constant 
is larger at the lower temperature. 



C. Fullerenes in Water 



For the purpose of modeling dilute solution of fullerenes in water, we have employed the following form for the 
solute-solvent interaction potential: [39] 



R-r 



R + r 



R-r 



R + r 



(26) 



where N=60, R=3.55 A, cr us =3.19 A, and e us =0.392 kJ/mol.Q] 
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FIG. 3. Methane-water radial distribution function at T=300 K and p=l g/cm 3 . Symbols are the simulation results for the 
methane- water radial distribution function with all-atom SPC/E water. [20I] and the line is from HNC theory for coarse-grained 
water described by the LJG potential. 



The fullerene-fullerene direct interaction potential is given by: |3£ 



m(r) = -a 



s(s-l) 3 s(s + l) 3 



s(s-l) 9 s(s + l) G 



2 



(27) 



where a=4.4775 kJ/mol, /3=0.0081 kJ/mol, and s = r/2R. 

In Fig . 6 we compare simulation results for the fullerene- water radial distribution function with all-atom SPC/E 
water, [20] and the HNC result for g(r) with coarse-grained LJG water. The agreement between the two distribution 
functions is satisfactory in general, although theory underestimates the heights of both primary and secondary peaks 
of g{r). Somewhat surprisingly, these discrepancies do not manifest themselves in the results for the solvent- mediated 
potential of mean force between two dilute fullerenes, as can be seen from Fig. 7, where we present simulation [l6j 
and theoretical results for AW(R). Theory and simulation are in good agreement throughout the entire range of 
separations studied, which could be a result of fortuitous cancellation of errors. 

In addition to the structural aspects of fullerene hydration, we have also studied dynamics of fullerene solutes 
in water. In particular, we have used the MCT approach outlined in the previous section to compute the diffusion 
coefficient of fullerene in water at thermodynamic state point T=300 K and p=l g/cm 3 ; we have obtained the value 
£>=1.73 nm 2 /ns. This can be compared to the simulation result[4(| of D=1.9 ± 0.4 nm 2 /ns, i.e. theoretical value is 
within the error-bars of the simulation result. 



IV. CONCLUSION 



In this paper we have performed a theoretical study of structural, thermodynamic, and dynamic properties of a 
coarse-grained water model based on a spherically symmetric LJG pair potential. For pure water, thermodynamically 
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FIG. 4. (a)PMF between two dilute methane solutes in coarse-grained LJG water; symbols are from MD simulations. [23l] and 
lines are from the integral equation theory, (b) Solvent-mediated PMF between two dilute methane solutes in coarse-grained 
LJG water; also shown are entropic and enthalpic terms. 



self-consistent HMSA integral equation theory was shown to give accurate results for microstructure and thermody- 
namics, while MCT theory was in good agreement with the simulation data for the self-diffusion coefficient, including 
its anomalous density dependence. 

In studying the structure of dilute solutions of hydrophobic solutes, anisotropic HNC theory was used to calculate 
the solute-solute PMF and the equilibrium constant for the conversion of the contact pair into the solvent-separated 
pair. Theory was shown to be in good agreement with the available simulation data for both methane-methane and 
fullcrcnc-fullcrcne PMF. It would be of interest to extend our theoretical approach to larger solutes, such as graphene 
plates; [13] this will be the subject of future research. 
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FIG. 6. Fullerene- water radial distribution function at T=300 K and p=l g/cm 3 . Symbols are the simulation results for the 
fullerene- water radial distribution function with all-atom SPC/E water. [2p|] and the line is from HNC theory for coarse-grained 
water described by the LJG potential. 
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FIG. 7. Solvent-mediated PMF between two dilute fullerene solutes in water. Symbols are the simulation results. [l6l] and lines 
are from integral equation theory. 



